% USRP Digtal Up Converter (DDC)  MatLab Simulation & Analysis

% By: Lee Kashka

% USRP DUC consists of 2 optional halfband filters with a 4 stage CIC 
% interpolating filter followed by AD9862 processing

function usrpduc(Fs_in)

%initialize everything for no Halfband filters
interp=32e6/Fs_in;
hb=0;
I_cic=interp;

%Determine if halfband filters should be used
if(rem(interp,4)==0)
    hb=2;
    I_cic=interp/4;
else
    if(rem(interp,2)==0)
        hb=1;
        I_cic=interp/2;
    end
end

% Number of CIC Sections
N_cic=4;

% HBF Interpolation Rate for each
I_hbf = 2;

%**************************************************************************
%************************* USRP HBF Filter Analysis************************
%**************************************************************************
% HBF Input Frequency is Fs_in MHz
Fs_hbf=Fs_in;

%coeff for large hb generated from round(2^18 * halfgen4(.7/4,8))/131072
Numerator = [-0.000816345214843750 0 0.00339508056640625 0 ...
             -0.00969696044921875 0 0.0225753784179688 0 ...
             -0.0465927124023438 0 0.0911941528320313 0 ...
             -0.188491821289063 0 0.628349304199219 1 ...
             0.628349304199219 0 -0.188491821289063 0 ...
             0.0911941528320313 ...
             0 -0.0465927124023438 0 0.0225753784179688 0 ...
             -0.00969696044921875 0 0.00339508056640625 ...
             0 -0.000816345214843750]; 
         
hbf = dfilt.dffir(Numerator);
info(hbf);
num  = get(hbf, 'Numerator');  % Get the numerator from the current filter.
hhbf=mfilt.firinterp(I_hbf,num);

%coeff for small hb generated from 2 * halfgen4(.75/8,2)
NumSmall = [-0.0815591197584173 0 0.578373182212901 ...
            1 0.578373182212901 0 -0.0815591197584173];
hbfs = dfilt.dffir(NumSmall);
info(hbfs);
num = get(hbfs, 'Numerator');  % Get the numerator from the current filter.
hhbfs=mfilt.firinterp(I_hbf,num);

%Gain Compensiation by dividing by 2 (corresponds to 1 bit shift)
hgain2=dfilt.scalar(1/2);

if(hb==2)
    hhbfnorm=cascade(hgain2,hhbf,hgain2,hhbfs);
    % Ploting Filter Frequency Response from 0 to Fs_cic/(4*R_cic) MHz
    h=fvtool(hhbfnorm,'Fs',Fs_hbf*4);
else
    if(hb==1)
       hhbfnorm=cascade(hgain2,hhbf);
       % Ploting Filter Frequency Response from 0 to Fs_cic/(2*R_cic) MHz
       h=fvtool(hhbfnorm,'Fs',Fs_hbf*2);
    end
end

%**************************************************************************
%*************************USRP CIC Filter Analysis*************************
%**************************************************************************

% CIC Output Frequency is 32 Mhz
Fs_cic=32e6;

hcic=mfilt.cicinterp(I_cic,1,N_cic);
% Print Filter Information
info(hcic);

% Ploting CIC Filter Frequency Response
%h=fvtool(hcic,'Fs',Fs_cic*I_cic);

% Do Gain Compensiation
hgain=dfilt.scalar(1/(I_cic*gain(hcic)));
hcicnorm=cascade(hgain,hcic);

% Ploting Filter Frequency Response After gain compensiation
h=fvtool(hcicnorm,'Fs',Fs_cic);

% Ploting Filter Frequency Response from 0 MHz to Fs_cic/2 MHz
%h=fvtool(hcicnorm,'Fs',Fs_cic);
axis([0 Fs_cic/(2*1e6) -120 5]);

%**************************************************************************
%************************* USRP AD9862 Filter Analysis*********************
%**************************************************************************

% Output Frequency is 128MHz
Fs_hAD=128e6;


%Calculated taps of 39 tap filter of AD9862 with 60 dB suppression
Numerator = [0.000031284734854580975 0.0022260936764527044 ...
   0.0034301041872895122 0.00079289970113994603 -0.0048353144417779652 ...
   -0.006111239917129716 0.0020408611126039246 0.011835550921197596 ...
   0.0079610077629940317 -0.010958704087977694 -0.02233902117979969 ...
   -0.0037725327677688883 0.030736045469094562 0.034209514286099453 ...
   -0.016820007868591794 -0.073149813195345761 -0.043746585382770882 ...
   0.10624491452824639 0.29473072539350925 0.38074483328498288 ...
   0.29473072539350925 0.10624491452824639 -0.043746585382770882 ...
   -0.073149813195345761 -0.016820007868591794 0.034209514286099453 ...
   0.030736045469094562 -0.0037725327677688883 -0.02233902117979969 ...
   -0.010958704087977694 0.0079610077629940317 0.011835550921197596 ...
   0.0020408611126039246 -0.006111239917129716 -0.0048353144417779652 ...
   0.00079289970113994603 0.0034301041872895122 0.0022260936764527044 ...
   0.000031284734854580975];     

hAD = dfilt.dffir(Numerator);
info(hAD);

num  = get(hAD, 'Numerator');  % Get numerator from the current filter.
hADf1=mfilt.firinterp(2,num);

%Calculated taps of 15 tap filter of AD9862 with 60 dB suppression
Numerator = [-0.0011444677289109315 -0.0014789503651362129 ...              
-0.0013415837550366906 0.000072851214910165318 ...            
 0.002909206612427241 0.0064808193320295008 ...              
 0.0091866691632470528 0.0089329163047706783 ...              
 0.00404597520318221 -0.0056112749164101069 ...              
-0.017815950400044765 -0.027959813617005874 ...               
-0.030085282626881701 -0.018816512544045975 ...               
 0.008424887484532885 0.049778949248305393 ...               
 0.098646698812048017 0.14505307407200793 ...                
 0.17836294702284217 0.19047408211595984 ...                
 0.17836294702284217 0.14505307407200793 ...                
 0.098646698812048017  0.049778949248305393 ...               
 0.008424887484532885 -0.018816512544045975 -0.030085282626881701 ...               
-0.027959813617005874 -0.017815950400044765 ...               
-0.0056112749164101069  0.00404597520318221 ...                
 0.0089329163047706783 0.0091866691632470528 ...              
 0.0064808193320295008 0.002909206612427241 ...               
 0.000072851214910165318 -0.0013415837550366906 ...              
-0.0014789503651362129 -0.0011444677289109315];  

hAD2 = dfilt.dffir(Numerator);
info(hAD2);

num  = get(hAD2, 'Numerator');
hADf2=mfilt.firinterp(2,num);

hADnorm=cascade(hADf1,hADf2);
% Ploting Filter Frequency Response from 0 to 128 MHz
h=fvtool(hADnorm,'Fs',Fs_hAD);

%**************************************************************************
%****************************Complete USRP DDC*****************************
%**************************************************************************


if(hb~=0)
    hduc=cascade(hhbfnorm, hcicnorm, hADnorm);
    h=fvtool(hhbfnorm, hcicnorm, hADnorm, hduc ,'Fs', ...
        [Fs_hbf*hb*2,Fs_cic,Fs_hAD, Fs_hAD]);
    axis([0 Fs_hAD/(2*1e6) -120 2]);
    legend(h,'HB','CIC','AD','HB+CIC+AD');
else
    hduc=cascade(hcicnorm, hADnorm);
    h=fvtool(hcicnorm, hADnorm, hduc ,'Fs',[Fs_cic,Fs_hAD, Fs_hAD]);
    axis([0 Fs_hAD/(2*1e6) -120 2]);
    legend(h,'CIC','AD','CIC+AD');
end

h=fvtool(hduc,'Fs',Fs_hAD);

end